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Abstract. The present paper pertains to corrections which are due to the presence 
of beam-limiting and beam-shaping devices in proton planning. Two types of 
corrections are considered; those originating from the nonzero thickness of such devices 
(geometrical effects) and those relating to the scattering of beam particles off their 
material. The application of these two types of corrections is greatly facilitated 
by decomposing their effects on the fluence into easily-calculable contributions. To 
expedite the derivation of the scattering corrections, a two-step process has been 
introduced into a commercial product which is widely used in planning (Treatment 
Planning System Eclipse*^™-'). The first step occurs at the beam-configuration phase 
and comprises the analysis of half-block fluence measurements and the extraction of 
the one parameter of the model which is used in the description of the beamline 
characteristics; subsequently, a number of Monte-Carlo runs yield the parameters of 
a convenient parameterisation of the relevant fluence contributions. The second step 
involves (at planning time) the reconstruction of the parameters (which are used in 
the description of the scattering contributions) via simple interpolations, performed on 
the results obtained during the beam-configuration phase. Given the lack of dedicated 
data, the validation of the method is currently based on the reproduction of the parts 
of the half-block fluence measurements (i.e., of the data used as input during the 
beam configuration) which had been removed from the database to suppress unwanted 
(block-scattering) contributions; it is convincingly demonstrated that the inclusion 
of the scattering effects leads to substantial improvement in the reproduction of the 
experimental data. The contributions from the block-thickness and block-scattering 
effects (to the fluence) are presented separately in the case of a simple water phantom; 
in this example, the maximal contribution of the two block-relating effects amounts to 
a few percent of the prescribed dose. 

PACS numbers: 87.55.Gh, 87.55.kh, 87.56.jk, 87.56.nk 
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1. Introduction 

To shape the beam so that it matches the characteristics of the specific treatment 
in radiation therapy (thus achieving the dehvery of the prescribed dose to the target 
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(tumour) and maximal protection of the surrounding healthy tissue and vital organs), 
beam-limiting and beam-shaping devices (BL/BSDs) are routinely used. Generally 
speaking, the beam is first restricted (in size) by the primary collimator, a beam- 
limiting device giving it a rectangular shape. The beam may subsequently encounter the 
multi-leaf collimator (MLC), which may be static or dynamic (i.e., undergoing software- 
controlled motion during the treatment session). More frequently than not, the desirable 
beam shaping is achieved by inserting a metallic piece (with the appropriate aperture 
and thickness) into the beamline, directly in front of the patient; this last beam-shaping 
device is called a patient collimator or simply a block. Being positioned close to the 
patient, the block achieves efficient fall-off of the dose (sharp penumbra) outside the 
target area. (The simultaneous use of MLC and block is not common.) 

To provide efficient attenuation of the beam outside the irradiated volume, the 
BL/BSDs are made of high-Z materials. The primary collimator and the MLC are fixed 
parts of the treatment machine, whereas the block not only depends on the particular 
patient, but also on the direction from which the target is irradiated. Therefore, there 
may be several blocks in one treatment plan, not necessarily corresponding to the same 
thickness or (perhaps) material. 

The presence of BL/BSDs in treatment plans induces three types of physical effects: 

a) Confinement of the beam to the area corresponding to full transmission (i.e., to 
the aperture of the device) . 

b) Effects associated with the nonzero thickness of the device (geometrical effects) . 

c) Effects relating to the scattering of the beam off the material of the device. 

Type- (a) effects (direct blocking of the beam) are dominant and have always been 
taken into account. The standard way to do this is by reducing the BL/BSD into 
a two-dimensional (2D) object (i.e., by disregarding its thickness) and assuming no 
transmission of the beam outside its aperture. Type-(b) and type-(c) effects induce 
corrections which, albeit at a few-percent level of the prescribed dose, may represent a 
sizable fraction of the local dose; due to their complexity and to time restrictions during 
the planning phase, these corrections have (so far) been omitted in clinical applications. 

The subject of the slit scattering in beam collimation was first addressed by Courant 
(1951). Courant extracted analytical solutions for the effective increase in the slit width 
(attributable to scattering) by solving the diffusion equation inside the collimator. To 
fulfill the boundary conditions, he introduced the negative-image technique, which was 
later criticised (e.g., by Burge and Smith (1962)). Despite the debatable usefulness of the 
practical use of Courant's work, that article set forth definitions which were employed 
in future research; for example, by categorising the scattered particles as: 

• those impinging upon the upstream face of the collimator and emerging from its 
inner surface (bore), 

• those entering the bore and scattering out of it, and 

• those entering the bore and leaving the downstream face of the collimator. 
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In the present paper, Courant's type-1 particles will correspond to our 'outer tracks' 
(OTs), type-2 particles to our 'bore-scattered inner tracks' (BSITs), and type-3 particles 
to our 'going-through inner tracks' (GTITs), see Fig. [H The tracks which do not hit 
the block will comprise the 'pristine' beam. 

Aiming at the determination of the optimal material for proton-beam coUimation 
and intending to provide data for experiments at the linear accelerator of the National 
Institute for Research in Nuclear Science (NIRNS) at Harwell, Burge and Smith (1962) 
re-addressed the slit-scattering problem and obtained a solution via the numerical 
integration of the diffusion equation. Burge and Smith reported considerable differences 
to Courant's results. In the last section of their article, the authors discussed alternative 
approaches to the slit-scattering problem. 

A Monte-Carlo (MC) method, as a means to study the coUimator-scattering effects, 
was introduced by Paganetti (1998). Using the GEANT code to simulate a proton 
beamline at the Hahn-Meitner Institut (HMI) in Berlin, Paganetti introduced simple 
parameterisations to account for the changes in the energy and angular divergence of the 
beam as it traverses the various beamline elements. Judging from Fig. 6 of that paper, 
one expects the scattering effects to be around the 1% level. In the last section of his 
article, Paganetti correctly predicted that 'Monte-Carlo methods will become important 
for providing proton phase-space distributions for input to treatment-planning routines, 
though the calculation of the target dose will still be done analytically.' Our strategy is 
similar to that of Paganetti: the dose, delivered by the pristine beam, will be corrected 
for scattering effects on the basis of results obtained via MC runs prior to planning 
(actually, when the particular proton-treatment machine is configured). 

Block scattering was also investigated, along the general lines of Paganetti's paper, 
in the work of van Luijk et al (2001). The (same version of the) GEANT code was used, 
to simulate a proton beamline at Kernfysisch Versneller Instituut (KVI) in Groningen, 
and the characteristics of the scattered protons were studied. To validate their approach, 
the authors obtained dose measurements for several field sizes and at several distances 
from the block. Unlike other works, the effect of scattering in air was also included 
in that study, and turned out to be more significant than previously thought (its 
contribution to the angular divergence of the beam exceeded 1 mrad per 1 m of air). 
One of the interesting conclusions of that paper was that the penumbra of the dose 
distribution is mostly accounted for by the lateral spread of the undisturbed beam; that 
conclusion somewhat allayed former fears that the extraction of the effective-source size 
(cr^) and of the effective source-axis distance (SAD), both obtained from measurements 
conducted in the presence of a block, might be seriously affected by block-scattering 
contributions. Finally, the paper addressed the asymmetries induced by a misalignment 
of the block, concluding that they might be sizable. 

In their recent work, Kimstrand et al (2008) put the emphasis on the scattering 
off the multi-leaf collimator (MLC). In their approach, they also obtain the values of 
the parameters (involved in their corrections) by using a MC GEANT-based method. 
In their Section 'Discussion and Conclusions', the authors, setting forth the future 
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perspectives, mention that 'a challenge ... is to implement collimator scatter for a pencil 
beam kernel dose calculation engine.' We agree on this being 'a challenge'. The authors 
then advance to pre-empt that '. . .the methods presented . . . can straightforwardly be 
applied to arbitrary shaped collimators of different materials, such as moulded patient- 
specific collimators used in passively scattered proton beams.' One should at least 
remark that their paper does not contain adequate information supporting the thesis 
that the proposed approach is of practical use in a domain where the execution time is 
a important factor. 

2. Method 

In the present paper, the standard ICRU (1987) coordinate system will be used; in 
beam's eye view, the x axis points to the right, y upward, and z toward the source. The 
origin of the coordinate system is the isocentre. 

2.1. Miniblocks 

Assuming vanishing transmission outside the aperture of the BL/BSD, all relevant 
contributions to the fluence involve 2D integrals over its area. The evaluation of such 
integrals is time-consuming; at a time when serious efforts are made to reduce the overall 
time allotted to each patient, the generous allocation of resources in the evaluation 
and application of corrections to the primary dose is unacceptable. To expedite the 
extraction of these corrections, one must find a fast way to decompose the effects of a 
2D object (e.g., of the aperture of the BL/BSD) into one-dimensional, easily-calculable, 
easily-applicable contributions. 

An example of one such 2D object is shown in Fig. [2j the area B represents the 
aperture (corresponding to full transmission of the beam) and is separated from the area 
A, representing the (high-Z) material, via the contour C. Full attenuation of the beam 
is assumed in the area A. The contour C and the 'outer' contour of the area A (which 
is not shown) may have any arbitrary shape (e.g., rectangular, circular, etc.); the only 
requirement is that the contour C be contained within the contour of the area A. 

Let us assume that the aim is to derive the influence of the BL/BSD at one point 
Q (Fig. [3]); the point Q is projected to the point P, on the BL/BSD plane. The point P 
may be contained in the interior or the exterior of the area B. 

2.1.1. The point P lies within the area B To enable the evaluation of the effects of the 
BL/BSD at the point Q, a set of directions on the BL/BSD plane, intersecting at 
the point P, are chosen; an example with = 4 is shown in Fig. |H The effects of the 
BL/BSD at the point Q will be evaluated from the elementary contributions of the line 
segments corresponding to the intersection of A^ straight lines with the contour C. These 
line segments, which are contained within the area B and are bound by the contour C, 
will be called miniblocks. The effect of the BL/BSD at the point Q will be evaluated by 
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averaging the elementary contributions of the miniblocks created around the point P. 
In Fig. m one of these miniblocks is represented by the line segment 5*1 5*2. Obviously, 
the number of miniblocks in each direction depends on the number of intersections of 
the straight line (drawn through the point P, parallel to the chosen direction) and the 
contour C. 

2.1.2. The point P lies outside the area B As in the previous section, a set of 
directions, intersecting at the point P, are chosen; an example with = 4 is shown in 
Fig. El Again, the effects of the BL/BSD at the point Q will be evaluated on the basis 
of the elementary contributions in these A^ directions. In Fig. [5l one of the miniblocks 
is denoted by 5*3 5*4. 

2.1.3. General remarks Evidently, the number of directions A^ is arbitrary. For fixed 
A^, the accuracy of the evaluation depends on the details of the contour C and on the 
proximity of the point P to it. The reliability of the estimation is expected to increase 
with A^. (The evaluation is exact for A^ ^ oo.) 

There is one difference between the cases described in Sections 12.1.11 and 12.1.21 In 
case that the point P lies within the area B, there will always be at least one miniblock 
per direction; if the point P lies outside the area B, there might be no intersections 
with the contour C in some directions (in which case, the corresponding elementary 
contributions vanish). 

Let us assume that the physical characteristics of the beam (lateral spread, angular 
divergence, energy, etc.) and the entire geometry (also involving the BL/BSD) are 
fixed. Apart from the original point Q, the elementary contributions will involve (the 
coordinates of) three additional points: the two end points of the particular miniblock 
and the point P. Due to the fact that the A^ directions are created around the point P, 
the two end points of every miniblock and the point P will always lie on one straight 
line. If the point P lies within the area B, it may lie within or outside a miniblock. 
On the other hand, if the point P lies outside the area B, it will always lie outside any 
corresponding miniblock. In any case, the elementary contribution of a miniblock will 
generally be a function of the distance of its two end points to the point P. Additionally, 
the elementary contribution of a miniblock to the point Q will also involve the distance 
of the point Q to the BL/BSD plane, denoted as z in Fig. [3l 

A few remarks are worth making. 

• The angle between consecutive directions may be constant or variable. If the angle 
is not constant, weights (to the elementary contributions) have to be applied. 

• The values of A^ in Sections l2.1.1l and l2.1.2l do not have to be the same. Furthermore, 
different values of A^ might be used for the points which are projected onto the 
interior (or exterior) of the area B, e.g., depending on the proximity of the point P 
to the contour C. 
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At present, the method of this paper apphes only to BL/BSDs with a constant 
aperture profile throughout their thickness; to derive the corrections in case of other 
shapes (e.g., for BL/BSDs with a tapered aperture profile), substantial modifications 
are required. The approach is applicable to any type of BL/BSD, be it the primary 
collimator, the MLC, or the block; to derive the corrections, the only input pieces 
are the physical characteristics of the material (of which the BL/BSD is made) and, 
naturally, geometrical details. From now on, however, we will restrict ourselves to the 
effects induced by the block, which (given its proximity to the patient) are expected 
to be of greater interest and importance in clinical applications. This choice, meant to 
emphasise importance, should not be seen as a restriction of the method. 



2.2. Block-thickness corrections 

A method for the evaluation of the thickness corrections was recently proposed by 
Slopsema and Kooy (2006); we will follow their terminology. Thus, 'downstream 
(upstream) projection' will indicate the projection of the downstream (upstream) face of 
the block on the {x,y) plane at the specified z position (depth). Similarly, the 'extension' 
of the block will correspond to the physical block translated in space (from its actual 
position to the specified depth, parallel to the central beam axis). 

In medical applications, the formal beam-shaper object (e.g., the so-called DICOM 
block), which is created at planning time to represent the physical block, comprises the 
projection of the downstream face (of the block) to a plane (perpendicular to the beam) 
at isocentre depth. The physical block (extension) is therefore obtained via a simple 
scaling involving the source-axis and block-isocentre (IBD) distances. Having retrieved 
the extension of the block, it is straightforward to obtain the downstream projection 
at any z; in order to obtain the upstream projections, one has to use, in addition 
to the aforementioned quantities, the block thickness (d). The relation between the 
DICOM block, the extension, and the downstream projection at a specified z position is 
shown in Fig. El The upstream projection is obtained by projecting the extension from 
depth z = IBD -|- d. The thin-block approximation, which is currently used in clinical 
applications when evaluating the dose, corresponds to d = mm (the downstream and 
upstream projections coincide at all z values). 

Assuming that the coordinates of the projected (downstream- or upstream-face, as 
the case might be) ends of a miniblock are denoted as xi and X2, the contribution of 
the miniblock to the fluence at a point Xp (lying on the straight line defined by Xi and 
X2) on the calculation plane is given by the formula 

F{x,) = \ I erf (^^) - erf (^-^) | , (1) 

where erf(x) denotes the error function and ai^2 stand for the rms (lateral) spreads of 
the beam at the specified depth. The quantities ai^2 are obtained via source mirroring 
according to the method of Slopsema and Kooy (2006); as points on the downstream 
or upstream faces of the block are used, the resulting a values in Eq. ([1]) are equal only 
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when one block face is involved in the mirroring process. 

One simple example of the projections on the calculation plane is shown in Fig. [71 
the block extension is contained within the downstream projection, which (in turn) is 
contained within the upstream projection. In reality, depending on the complexity of 
the shape of the block aperture and on the relative position of the central beam axis, 
these three contours might intersect one another. The central beam axis intersects the 
calculation plane at the origin of the {x,y) coordinate system. One miniblock is shown 
(the line segment contained within the block extension), along with two points, one 
within the miniblock (P), another outside the miniblock (P'). The fluence contributions 
to both points may be evaluated by using Eq. ([T]) with the appropriate xi^2 and ai^2 
values. 

The essentials for the evaluation of the contribution of a miniblock to the fluence 
at a specified point on the calculation plane are to be found in Fig. [HI Although, in 
the general case, the miniblock does not contain the intersection of the x and y axes of 
Fig. [71 all lengths (which are important to our purpose) scale by the same factor, thus 
enabling the simplified picture of Fig. [81 

According to Slopsema and Kooy (2006), the source is mirrored onto the calculation 
plane by using points either on the downstream or on the upstream face of the block. 
The values of the mirrored-source size, corresponding to these two options, are given by 

IBD - z 
= SAD -IBD 

and 

IBD-z + d 

''""''^SAD-IBD-rf- ^"^^ 
The coordinates of the projected points Ui, Di, D2, and U2 are obtained on the 
basis of Fig. [HI via simple operations. The last item needed to derive the contributions 
to the fluence may be found on page 5444 of the paper of Slopsema and Kooy (2006): 
'Protons whose tracks project inside the aperture extension onto the plane of interest 
see the upstream face as the limiting boundary. Only for protons whose tracks end 
up outside the aperture on the plane of interest is the downstream face the limiting 
aperture boundary.' In fact, this statement applies to the case of a half-block (one-sided 
block). The modification, however, in case of a miniblock is straightforward. 

• Points within the extension of the aperture see the upstream face of the miniblock 
as the limiting boundary. 

• Points outside the extension of the miniblock see one upstream and one downstream 
edge as limiting boundaries. 

Obviously, in order to evaluate the fluence at the point P (Fig. [8|), one has to use the 
coordinates of the projected points Ui and U2, along with (j„ in Eq. ([T]). On the other 
hand, to evaluate the fluence at P', one has to use Ui along with cr„ (contribution of 
the 'left' part of the miniblock) and D2 along with ad (contribution of the 'right' part 
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of the miniblock). This simphfied picture, featuring what a point 'perceives' as hmiting 
boundaries, suffices in obtaining the appropriate fluence contributions. 

Finally, the number of contributions which a point will receive depends on the 
geometry and on the shape of the block. One example of a point receiving two 
contributions in a given direction is shown in Fig. [9l the point P lies within the extension 
of the miniblock on the left and outside the extension of the miniblock on the right. 
Additionally, one might have to deal with a block with more than one apertures (hence 
contours). To determine in the present work the appropriate thickness corrections 
for arbitrary block-aperture shape and beamline geometry, dedicated software was 
developed. 

2.3. Block-scattering corrections 

Before entering the details of the derivation of the scattering corrections, it is worth 
providing a concise outline of our approach. First of all, our aim is to obtain a fast and 
reliable solution to the scattering problem; exact analytical solutions are welcome as 
long as they fulfill this requirement. Second, the solution has to be general enough for 
direct application to all proton-treatment machines. 

To expedite the derivation and application of the corrections at planning time, we 
will introduce a two-step approach. 

• All the parameters which are independent of the specific details of plans will be 
evaluated during the beam-configuration phase; given a fixed (hardware) setup, 
every treatment machine is configured only once. 

• The corrections for each particular plan will be derived (at planning time) from 
the existing results (i.e., those obtained at beam-configuration phase) via simple 
interpolations. 

Given that the physics of multiple scattering is known, it is possible to obtain the 
exact solution for the relevant beam properties (lateral spread, angular divergence) from 
the beamline characteristics of each treatment machine. However, it is unrealistic to 
introduce a dedicated process for each supported machine, especially in order to derive 
corrections to the delivered dose. Additionally, if a dedicated (per-case) approach is 
implemented in a software product which is intended to support a variety of machine 
manufacturers, one has to be prepared to allot the necessary resources whenever a new 
product appears. To avoid these problems and to retain the generality of the approach, 
one has no other choice but to introduce a simple, adjustable model to account for the 
beam optics. The model of the present paper has only one parameter, which will be 
fixed from half-block fluence measurements. 

The parameters which achieve the description of the various distributions of the 
scattered protons are determined via MC runs. These runs take account of the variability 
in the block material, block thickness, incident energy (energy at nozzle entrance), and 
nozzle-equivalent thickness (NeT) for all the options (combinations of the hardware 
components of the beamline, leading to ranges of available energies and of NeTs, as well 
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as imposing restrictions on the field size) for which a machine is configured. To enable 
the easy use of the results, the output is put in the form of expansion parameters in two 
geometrical quantities which are involved in the description of the scattering effects. 

The scattering corrections for all the blocks in a plan are determined (at planning 
time) from the aforementioned results via simple interpolations. The application of the 
corrections involves the concept of miniblocks, as they have been introduced in Section 
12. II of this paper. 



2.3.1. Modelling of the beam One model which is frequently used in beam optics 
features the bivariate Gaussian distribution in the lateral direction y (distance to the 
central beam axis) and the (small) angle 9 (with respect to that axis). (Rotational 
symmetry is assumed here.) 

.... 1 / Ae'-2Bye + Cy\ 

^^^'^^ = D ) ' 

with 

D = AC-B'^. (5) 

The parameters A and C represent twice the variance in y and 6*, respectively. The {yfi) 
correlation is defined as 

The quantity p (which is bound between —1 and 1) is a measure of the focusing in 
the beam. Positive p values indicate a defocusing system, negative a focusing one; this 
becomes obvious after one puts Eq. (jlj) in the form 

1 fl2 1 ill — o—9\^ 

/(y,g) = ^^exp(-f^) .. ^^P(- 9 2n J - 



We now touch on the variation of A, B, and C along the beam-propagation 
direction. Assuming that the quantities Aq, Bq, and Cq denote the corresponding values 
at isocentre depth and that the beam propagates in air (without scattering). A, B, and 
C at distance z from the isocentre (see Fig. [6]) are given by the expressions 

A{z) = Ao- 2Boz + Coz^ , B{z) = Bo - Cqz , and C{z) = Cq . (8) 

With these transformations, the joint probability distribution of Eq. (jlj) is invariant 
under translations in z. In case that the beam traverses some material, Eqs. (IE]) have 
to be modified accordingly, to take account of the beam broadening due to multiple 
scattering. 



2.3.2. Simplified parameterisation of the heamline The accurate modelling of the beam 
may be obtained on the basis of formulae such as those given in the previous section. 
In Section 12.31 however, we reasoned that a simplified parameterisation of the beamline 
is desirable; one additional argument may be put forth. 
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Currently, as far as proton therapy is concerned, four dose-delivery techniques 
are in use: single-scattering, double-scattering, uniform-scanning (formerly known as 
wobbling), and modulated-scanning (formerly simply known as scanning). In the 
modulated-scanning technique, magnets deflect a narrow beam onto a sequence of 
pre-established points (spots) on the patient (for pre-determined optimal times), thus 
'scanning' the (cross section of the) region of interest. Uniform scanning involves the 
spread-out of the beam using fast magnetic switching. The broadening of the beam in 
the single-scattering technique is achieved by one scatterer, made of a high-Z material 
and placed close to the entrance of the nozzle. Currently, the most 'popular' technique 
involves a double-scattering system. 

In a double-scattering system, a second scatterer is placed downstream of the first 
scatterer in order to achieve efficient broadening of the beam; studies of the effects of 
the second scatterer may be found in the literature, e.g., see Takada (2002) and more 
recent reports by the same author and Gottschalk. The second scatterer is usually 
made of two materials: a high-Z (such as lead) material at the centre (i.e., close to the 
central beam axis), surrounded by a \ow-Z (such as aluminium, lexan, etc.) material 
(which is frequently, but not necessarily, shaped as a concentric ring) . The arrangement 
produces more scattering at the centre than the periphery, leading (after sophisticated 
fine-tuning) to the creation of a broad flat field at isocentre. 

To simulate the effect of the second scatterer in the present work, {y,6) events 
arc generated at z = SAD as follows. The variable y is sampled from the Gaussian 
distribution with mean and variance 0"^ (the source-size calibration must precede this 
step); as depends on the incident energy and NeT. To account for the lateral limits of the 
beam, y is restricted within the interval [—Rl,Rl]: where the characteristic length Rl 
is taken herein to be the radius of the second scatterer. The variable 9 is first sampled 
from the Gaussian distribution with mean and a y-dependent variance according 
to the formula 

a,(t) = (7,(l)((l-A) |t|+A), (9) 

where t denotes the lateral position as a fraction of Rl, \ t | = | y/ Rl |< 1- To account for 
the y-dependent bias in 9, we then use the transformation 9 — > 9 -\-isin~^ {y / L) ~ 9+y/ L, 
where L stands for the distance between the first scatterer and the source. Obviously, 
A denotes the ratio of two ae values, i.e., the value at the centre of the second scatterer 
over the one at the rim; A is the only free parameter of the model. The value of (7^(1) 
is obtained from the incident energy and NeT; the angular divergence of the beam at 
nozzle entrance is (currently) assumed to be 0. 

It is not difficult to prove that, given the aforementioned rules for generating the 
{y,9) events, the open-field fluence at a lateral position y at depth z < SAD is given by 

^^y^ = 2nasiSAD-z) L ^ ^"K-^)-P( ) ' (1°) 

where to = 1. It has to be emphasised that this definition of the fluence does not 
involve an overall ^-factor (r being the distance between the calculation plane and 
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the source); therefore, this formula is compatible with the format in which the lateral 
fluence measurements, used during the beam-configuration phase, appear. 

It can be shown that the only modification in case that a half-block is inserted 
into the beamline (e.g., as shown in Fig. [1]) involves the upper limit to of integration in 
Eq. (ITOl) : instead of to = 1, one must now use 

. f 6(SAD - z) - y(SAD - IBD) 1 

if to > ~1) or otherwise to = —1. 

Although the method of the present paper was originally developed for the double- 
scattering technique, it is also applicable in single scattering and uniform scanning; in 
both cases, one simply has to fix the parameter A to 1. As the method is applicable 
only in case of broad fields, it has no bearing on modulated scanning. 

2. 3. 3. Multiple scattering through small angles The elements needed for the description 
of the passage of particles through matter may be found (in a concise form) in Yao et 
al (2006), starting on page 258; most of the deflection of a charged particle traversing a 
medium is due to Coulomb scattering off nuclei. Despite its incompleteness (e.g., see the 
discussion in the GEANT4 physics-reference manual, section on 'Multiple Scattering', 
starting on page 71), the multiple-scattering model of Mohere is used here. The large- 
angle scattering is not taken into account; the angular distribution of the traversing 
beam is assumed Gaussian. 

Highland's logarithmic term, appearing in the expression of Oq, will be approximated 
by a constant factor involving the block thickness; a similar strategy was followed 
in Gottschalk et al (1993). The Lynch-and-Dahl (1991) values will be used in the 
expression for ^o^ 

»„(,) = lMM£l/Z(i + o.0381n(A)), (12) 

pep V ^0 ^ ^0 ^ 

where q denotes the depth along the original trajectory, (3c the velocity and p the 
momentum of the proton, and Xq the radiation length in the material of the block (for 
a convenient parameterisation of Xq, see Yao et al (2006), page 263). 

Equation (fT2l) applies to 'thin' targets. For 'thick' targets, the dependence of the 
proton momentum on the depth q must be taken into account. Omitting the logarithmic 
term on the right-hand side, one may put Eq. f[T2|) in the form 

0oiq) = 2fiq)^, (13) 

where 

, 13.6MeV/l /"? dt Ni/2 



On the practical side, Eq. (fT3l) with a constant (or at least 'not too complicated') /(g) 
factor would be attractive as one would then be able to obtain fast analytical solutions 
for the propagation of a simulated track inside the material of the block. Therefore, 
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it is worth determining the range of thickness values within which the constancy of 
the f{q) factor remains a reasonable assumption. The direct comparison with the data 
of Gottschalk et al (1993) revealed that the 'thick-target' corrections are unfortunately 
indispensable at depths exceeding about 50% of the range of 158.6-MeV protons, incident 
on a variety of materials. 

To abide by the original goal of obtaining a fast solution, we had to follow an 
alternative approach (to that of using Eq. f|T4|) ). by parameterising the g-dependence of 
the /(g) factor in a simple manner; at present, the best choice seems to be to make use 
of the empirical formula 



~ 1.75/J2 

where R denotes the (energy-dependent) range of the incident proton in the material of 
the block and 



/ is now a constant, depending only on the initial value of (3. The validation of Eqs. (ITSll 
and (JT6l) was made on the basis of a comparison with the experimental data of Gottschalk 
et al (1993), namely the measured 6m values of that paper. Good agreement with the 
data was obtained for four materials which are of interest in the context of the present 
study (carbon, lexan, aluminum, and lead). Only at one entry (one of the largest depth 
values in lead, i.e., the measurement at q/R = 0.97548), was a significant difference (of 
slightly less than 20%) found; the origin of that difference was not sought. 

2.3.4- Details on the generation of the MC events Figure fTOl shows an example of a 
trajectory of one proton inside a block. The incident proton, an OT in this figure, 
hits the upstream face of the block at angle 9 with respect to the beam axis. Two 
new in-plane variables are introduced to describe the kinematics at depth q (along the 
direction of the original trajectory): the deflection r (off the original-trajectory course) 
and the angle (with respect to the direction of the original trajectory). Although the 
proton moves in an irregular path inside the block, the 'history' of the actual motion 
will be replaced by a smooth movement leading to the same value of r at q. This 
'smooth-deflection' approximation will enable the association of the energy loss in the 
material of the block with the doublet of (g,r) values. Since the path length, which 
is calculated in this approximation, is an underestimate of the actual path length, a 
constant [§| conversion factor (true-path correction) of 0.9 has been used; there is some 
arbitrariness concerning this choice, yet it appears to be reasonable. It needs to be 
stressed that, in Fig. [10], x denotes the direction of the beam propagation; the auxiliary 
coordinate system introduced in this figure should not be confused with the coordinate 
system of Fig. [6l which is the formal one in medical applications. 



^o(g) 



2/v/g 



(15) 



1 - 




(16) 



§ In reality, the true-path correction is expected to be energy-dependent. 
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In the generation of the MC events, the suggestion of Yao et al (2006), page 262, for 
the quantities r and (p has been followed. Two independent Gaussian random variables 
(^1,2:2) with mean and variance 1 are first created in each event (track hitting the 
block). The quantities r and are expressed in terms of (-21,-22) as follows: 

r(g) = (^ + ^2j^- (17) 

and 

0(g) = z^eoiq) , (18) 

where OqIq) is taken from Eq. (fT5|) . 

The values of the doublet (-21,-22) fix the dependence of r and on g in each generated 
event. It may then be determined (either analytically or numerically) whether the 
particular track leaves the block. Finally, for those tracks leaving the block, a simple 
rotation yields the coordinates of the exit point in the {x,y) coordinate system of Fig.fTOl 
The energy of the leaving proton is determined from its residual range (original range 
minus the actual path length inside the material of the block). 

2. 3. 5. The lateral fluence distributions of the scattered protons The definitions of the 
three types of scattered protons have been given in Section [U in Fig. [H one obtains a 
rough schematic view of the corresponding contributions to the fluence. In this section, 
we will introduce convenient forms to parameterise the lateral fluence distributions of 
the scattered protons. 

As far as these distributions are concerned, the uninteresting offset h (lateral 
displacement of the block) in Fig. [1] will be omitted. Therefore, for the needs of this 
section, ?/ = mm at the extension of the block in Fig. [H that is, not at the position 
where the central beam axis intersects the calculation plane. To obtain the lateral 
fluence distributions, literally corresponding to Fig. [1], one has to take the offset h into 
account. In the formulae below, the placement of the half-block is assumed as shown in 
Fig. [1] (i.e., extending to +00). 

The empirical formula for the description of the lateral fluence distribution of the 
OTs reads as 

f{y) = aeM-y/P){y + l)\ (19) 

where y < 0. Four conditions for the parameters a, /3, and 7 must be fulfilled: a > 0, 
/9 < 0, 7 < 0, and 2/3 - 7 < 0. 

The same empirical formula is used in the parameterisation of the lateral fluence 
distribution of the BSITs; the four aforementioned conditions also apply. (The resulting 
optimal parameter values are, of course, different.) 

The optimal description of the (broader distribution of the) GTITs is achieved on 
the basis of a Lorentzian 
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multiplied by the asymmetry factor 

2 



l + exp(-2^)' 



to account for the observed skewness of the lateral fluence distribution (toward positive 
y values). Three conditions must be fulfilled: a > 0, /? > 0, and 7 > 0. 

The lateral fluence distribution of the pristine tracks is fitted by using the standard 
formula 



where a > and 7 > 0; ^(x) is the cumulative Gaussian distribution function. Only 
the factor a will be retained from the fits to the pristine-beam data, to be used in 
the normalisation of the fluence corresponding to each of the three types of scattered 
protons. Expressing the contributions of the scattered protons as fractions of the 
pristine-beam fluence enables the efficient application of the corrections at planning 
time. 

From all the above, it is evident that the description of the contribution to the 
fluence (at fixed z) of any of the three types of the scattered protons is achieved on the 
basis of three parameters; hence, there are nine parameters in total. At the end of each 
cycle (comprising a set of MC runs for a number of b values, for a given energy-NeT 
combination), each of these nine parameters Pi is expanded in terms of b and z, using 
the quadratic model 



The final results are the 56 coefficients aij, i G {1,...,9}, j G {1,...,6}, obtained 
at several energy-NeT combinations for the option of the machine which is under 
calibration. At planning time, the values of the parameters Pi are reconstructed from 
the existing results (via simple 2D interpolations in the incident energy and NeT), 
the values of b (corresponding to the particular miniblock which is processed), and z 
(corresponding to the calculation plane which is processed). 

2. 3. 6. A summary of the approach Let us assume that one option of a selected machine 
has been chosen for calibration (at beam-configuration phase). All half-block data 
(which is contained in that option) is used in the determination of the parameter A on 
the basis of an optimisation scheme (featuring the C++ implementation of the standard 
optimisation package MINUIT of the CERN library), along with Eq. ffTOl) with the 
value of Eq. ^ % 

Representative incident-energy and NeT values are chosen from the ranges of values 
associated with the option which is under calibration. For each acceptable energy-NeT 

II For the determination of the parameter A, one could also use the measurements of the open-field 
fluence, along with Eq. pO|l with io = 1- However, what is habitually called 'open-field measurements' 
in the field of radiation therapy corresponds to beams which have already been restricted in size by the 
primary collimator. 




(21) 



Pi = an + aj2& + CLi^z + ai^bz + 0^56^ + OiQz'^ . 



(22) 
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combination, a number of MC runs are performed, each corresponding to one b value. 
In each of these MC runs, events (i.e., {y,0) pairs, each corresponding to one proton 
track) are generated according to the formahsm developed in Section I2.3.2[ The value 
of the parameter A, obtained from the half-block data (for the option in question) at 
the previous step, comprises (i.e., apart from geometrical characteristics of the machine 
which is configured) the only input to these MC runs. The resulting tracks are followed 
until they either hit the block or pass through it. The tracking of the protons inside the 
material of the block is done according to Section [2331 finally, these tracks either vanish 
(being absorbed in the material of the block) or emerge from it (bore, downstream face) 
and deliver dose. 

The tracks which emerge from the block are properly flagged (OTs, BSlTs, or 
GTITs) and their contributions to the fluence at a number of z depths are stored 
(histogrammed). Fits to these distributions, using the empirical formulae of Section 
I2.3.5[ lead to the extraction of the parameters achieving the optimal description of the 
stored data. After the completion of all the runs for all the chosen values of b, the entire 
set of the parameter values is subjected to fits by using Eq. fj22|) for each parameter 
separately. Finally, the coefficients aij, i G {1,...,9}, j G {1,...,6}, appearing in 
Eq. (122|) . are stored in files along with the values of the incident energy and NeT 
corresponding to the particular MC run; these output files will comprise the only input 
when the corrections to a particular plan will be derived (planning time). 

The procedure above is repeated until all options of the given proton-treatment 
machine have been calibrated. The variability in the material and in the thickness of 
the block is also taken into account in the current implementation (by looping over 
those combinations requested by the user0). Finally, it is worth repeating that this 
time-consuming part of obtaining the files containing the a^- coefficients (a question 
of a few hours per option) has to be performed only once, when the proton-treatment 
machine is configured. 

At planning time, the pristine-beam fluence is calculated first. The beam-scattering 
corrections are then obtained (i.e., if they have been requested) after employing a number 
of elements developed in the course of the present section: 

• The concept of miniblocks introduced in Section 12. 1[ 

• The reconstruction of each of the nine parameters, used in the description of the 
scattered protons, at a few z values, on the basis of Eq. ( l22l) from the results 
pertaining to the option selected in the treatment plan. The appropriate file 
(corresponding to the block material and thickness in the plan) is used as input. The 
final results for the various parameters are obtained via simple 2D interpolations 
in the incident energy and NeT. 

• The empirical formulae of Section I2.3.5[ 

^ At present, five materials are supported: brass, cerrobend, nickel, copper, and lead; this list is easily 
expandable. The block-thickness values are taken from the option properties of the machine which is 
configured. 
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• A simple (linear) interpolation to obtain the corrections at all depths z in the plan. 

The final step involves the application of the corrections to the pristine-beam fluence. 

The break-up of the task of determining the scattering corrections into two steps, as 
described in this section of the paper, enables the minimisation of the time consumption 
during the evaluation and application of these corrections in proton planning; of central 
importance, in this respect, is the concept of miniblocks. 

3. Results 

3.1. Machine configuration 

The measurements, which are analysed in this paper, have been obtained at the Proton 
Therapy Centre of the National Cancer Centre (NCC), South Korea. The first report on 
the clinical commissioning and the quality assurance for proton treatment at the NCC 
appeared shghtly more than one year ago, see Lee (2007). 

The NCC proton-treatment machine has been manufactured by Ion Beam 
Applications (IBA), Louvain-la-Neuve, Belgium. Its nominal SAD is 2300 mm and the 
distance between the first scatterer and the isocentre is 2792 mm. The double-scattering 
technique currently supports eight options with incident energies ranging from 155 to 
230 MeV. 

3.1.1. Half-block fluence measurements All half-block measurements have been taken 
in air, using IBD = 250 mm. The lateral displacement of the 65mm-thick brass block, 
used in these measurements, was b = —50 mm; the block was positioned opposite to 
what is shown in Fig. [T], i.e., blocking to beam from —50 mm to (theoretically) —oo. 
To apply Eq. ( JTOl) with the to value of Eq. (ITTl) . the y axis was inverted, as a result 
of which the b value of -|-50 mm was finally used in Eqs. (fTOj) and ffTTj) . Each option 
of the double-scattering technique at the NCC comprises 15 energy-NeT combinations; 
in each of these combinations, the lateral fluence distributions have been obtained at 
four z positions, namely at z = —300, —150, 0, and 150 mm. One example of these 
profiles is given in Fig. [TTl It is worth mentioning that each profile had been separately 
normalised (during the data taking) by setting the corresponding average value of the 
fluence 'close' to the central beam axis to (the arbitrary value of) 100 0; unfortunately, 
the individual normalisation factors are not available. The 'ears' of the distribution for 
the data set at z = 150 mm, which are presumably due to block scattering, have been 

+ At present, the analysis of the half-block fluence measurements is tedious. Unfortunately, there is no 
standard format in which these measurements, taken at the various treatment centres and machines, 
have to appear; in fact, there is complete freedom (during machine configuration) in the choice of the 
distance between the downstream face of the block and the isocentre, in the thickness and the lateral 
displacement of the block, and in the normalisation factors used in the output files. A more serious 
drawback is that the measurements are frequently (luckily, not in the NCC case) shifted laterally, so 
that the 50% of the corresponding maximal fluence (of each set, separately) be brought to y = mm; 
thus, an important offset is irretrievably lost. 
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removed via software cuts. To avoid fitting the noise, fluence values below 10% of the 
maximal value in each data set were not processed. The block-thickness effects were 
removed from the data prior to processing. 

Before advancing to the results, one important remark is prompt. A significant 
deterioration of the description of the data in the last of the options (option 8) was 
discovered during the analysis; at present, it is not entirely clear what causes this 
problem. It seems that our beamline model leads to a systematic overestimation of 
the penumbra in option 8. It has to be mentioned that the data in that option yields an 
unusually small spot size (a^), about 3 times smaller than the typical values extracted in 
the other seven options. Option 8 is the only one in which the IBA second scatterer SS3 
(which is admittedly of peculiar design 0) is used. In all probability, the problems seen 
in option 8 originate from the shape of the second scatterer, namely from the fact that 
the energy loss in its material is not kept constant radially (actually, it is 'discontinuous' 
at t ~ 0.9). In order to investigate the sensitivity of our conclusions to the treatment of 
option 8, we will perform the analysis both including in and excluding from the database 
the measurements of that option. 

3.1.2. Extraction of parameter A The extracted values of the parameter A for all the 
energy-NeT combinations of all double-scattering options of the NCC machine are shown 
in Fig. [T2I The uncertainties in the case of option 8 are (on average) larger than in the 
other options. The variability of the values within each option is due to the fact that, as 
a result of the numerous assumptions made to simplify the problem, A (which should, 
in principle, characterise only the second scatterer) was finally turned into an effective 
parameter. To decrease the 'noise' seen in Fig. [121 the weighted average of the extracted 
A values within each option was finally used in the ensuing MC runs for all energy-NeT 
combinations in that option (see Table [1]). 

3.1.3. Monte-Carlo runs 50 million events have been generated in each energy-NeT 
combination, at each of three h values (20, 40, and 60 mm). The lateral fluence 
distributions have been obtained at 17 positions in depth, from z = 100 to —100 mm. 
The parameters of these distributions for the three types of scattered protons have been 
extracted using the formulae of Section I2.3.5[ In all cases (i.e., including option 8), 
the description of the data was good; the reduced values (x^/NDF, NDF being the 
number of degrees of freedom in the fit) came out reasonably close to 1. Figures [13] and 
[TH show typical lateral fluence distributions for outer and inner tracks, respectively. 

* The amount of lead, used at the centre of the SS3, is smaller than in the cases of the SSs and SS2 
scatterers which are used in options 1-7. Furthermore, the thickness of lexan abruptly increases close 
to the rim of this scatterer. Evidently, the amount of material used in the SS3 cannot provide efficient 
broadening of the beam; adding material would imply smaller beam energy at nozzle exit, at a time 
when the emphasis in this option is obviously put on the high end of the energy spectrum. To be 
compatible with the requirement for flatness of the resulting field at isocentre, a smaller maximal field 
diameter (140 versus 220 to 240 mm of the other options) in the clinical application of option 8 had to 
be imposed. 
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Figures [TMT71 obtained with a lateral block displacement of 6 = 20 mm, show the 
scattering corrections (to be applied to the lateral fluence distributions of the pristine 
tracks) at three z positions around the isocentre (100, 0, and —100 mm, respectively). 
As expected, the distributions broaden when receding from the block; the mode of the 
fluence contribution of the OTs moves about 1 mm away from the block extension for 
every 10 mm of depth, thus indicating an 'average' exit angle (to the bore) of about 6°. 
Concerning their magnitude, the corrections generally amount to a few percent of the 
corresponding pristine fluence for the typical distances involved in clinical applications. 

It is now time to enter the subject of the energy loss of the scattered protons. To 
a good approximation, one may assume that the energy distributions of the scattered 
protons depend only on the ratio uj = E/E^g^^, where E is the energy of the scattered 
proton and ^'max denotes the energy of the pristine beam (i.e., the energy at nozzle 
exit). In their study, Kimstrand et al (2008) made the same observation. The energy 
distributions of the scattered protons were investigated in the case of three energy-NeT 
combinations of the NCC machine: (158.42 MeV, 125.3 mm), (182.99 MeV, 65.6 mm), 
and (229.75 MeV, 40.7 mm). These combinations have been selected from the data sets 
of options 1, 5, and 8, respectively, in which different second scatterers are employed; 
the corresponding energies -Emax are: 79.5, 150.2, and 212.8 MeV. The results have 
been obtained using a 65-mm thick brass block. Figures [18] and [19] show the energy 
distributions of the scattered protons as functions of the ratio oj; instead of referring 
explicitly to each energy-NeT combination, we will simply use the option number (or 
-Emax) to identify the results. The following remarks are worth making. 

• At low and moderate values of the exit energy £^niax; the energy distributions of 
the two types of inner tracks are similar, following the f{oj) = 2uj probability 
distribution. At the high end of the energies used in the NCC machine, the energy 
distribution of the BSITs departs from this simplified picture, attaining a peak 
close io UJ = 1; this is due to BSITs which almost 'brush' the surface of the bore. 
The energy distribution of the GTITs remains unchanged. Nevertheless, to retain 
simplicity, we will assume that the energy distribution of all inner tracks follows 
the formula f{uj) = 2uj. 

• The energy distribution of the OTs is strongly smeared toward low uj values. The 
OTs lose a significant amount of their energy when traversing the material of the 
block; as shown in Fig. [18], their energy distributions peak around uj ~ 0.2 to 0.3. 
To fit the energy distribution of the OTs, we used the empirical formula 

f{uj) ~ A/^exp(-?7u;) , 

where the parameter rj turns out to be around 2 to 2.5; the constant of 
proportionality is obtained from the normalisation of the probability distribution. 

Courant's effective-size corrections, corresponding to the three aforementioned 
cases, are: 0.31, 0.87, and 1.57 mm. At small aperture sizes, the dominant contribution 
to the fluence (of the scattered protons) originates from OTs; the inner tracks dominate 
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at high aperture values. The crossover is energy-dependent, ranging from about 15 
(option-1 result) to 47 mm (option-8 result). 

The absolute yields of the different types of the scattered protons are given in Table 
[2] for an aperture size (i.e., the diameter, assuming circular shape) of 100 mm. We 
observe that the absolute yield of the OTs increases by a factor of about 2 between 79.5 
and 150.2 MeV, and by more than 4 between 150.2 and 212.8 MeV. At 79.5 MeV, the 
BSITs account for more than 50% of the total yield of the scattered protons. However, 
the importance of the BSITs diminishes with increasing energy, reaching the level of 1 
out of 4 OTs at 212.8 MeV. On the contrary, the yield of the GTITs flattens out at 
about 0.75 per OT. One might conclude that the BSITs seem to be more important at 
low energies and the OTs at high energies. The ratio of the yields of the GTITs and 
OTs is almost constant, varying only from 2/3 to 3/4 as the energy increases from 79.5 
to 212.8 MeV. 

3.2. Verification of the method 

The verification of the method of the present paper should obviously involve the 
reproduction of dedicated dose measurements obtained in some material which, as far 
as the stopping power for protons is concerned, resembles human tissue, e.g., in water. 
The measurements should cover the region around the isocentre, where the tumour is 
usually placed, and, in order that the approach be validated also in the entrance region, 
they should extend to small distances from (the downstream face of) the block. Finally, 
the method must be validated for a range of depth values (associated with the energy 
at nozzle exit). 

At present, given the lack of dedicated dose measurements, the only possibility 
for verification rests on re-using the calibration data (i.e., the half-block fluence 
measurements described in Section 13.1.11) . We are aware of the fact that using the 
same data for configuring a system and for validating its output does not constitute 
an acceptable practice. However, since parts of the data (i.e., the areas which are 
obviously contaminated by the block-scattering effects) had been removed from the 
database before extracting the A values, this approach becomes a valid option. Luckily, 
as far as the validation of the scattering corrections to the fluence is concerned, our 
interest lies in (the reproduction of) those excluded areas; naturally, for the rest of the 
measurements (i.e., for those which were used in the determination of the A values), 
it has to be verified that the quality of the description of the experimental data is not 
impaired by the inclusion of the block-scattering corrections. 

A typical reproduction of the measurements is given in Figs. [20] and [211 the data 
correspond to the first energy-NeT combination of option 1 of the NCC machine, taken 
at z = 150 mm (i.e., 100 mm away from the block). Shown in Fig. [20] are the lateral 
fluence measurements (continuous line) along with the MC data corresponding only to 
the pristine beam; the effects of the scattered protons are added to the pristine-beam 
fluence (resulting in what will be called henceforth 'total fluence') in Fig. [21] On the 
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basis of the visual inspection of these two figures, there is no doubt that the quahty 
of the reproduction of the measurements in the latter case (i.e., when including the 
block-scattering effects) is superior. 

We will next investigate the goodness of the reproduction of all measurements on 
the basis of a commonly-used statistical measure, e.g., of the standard function. 
Alternative options have been established (e.g., the 7- index approach of Low et al 
(1998)), but have not been tried in this work. One has to bear in mind that the 
block-scattering contributions are larger at small distances from the block and in the 
area neighbouring its extension; at large distances, the distributions of the scattered 
protons broaden as a result of the angular divergence of the scattered beam and, less 
importantly, of scattering in air (an effect which has not been included in this paper). 
Evidently, the assessment of the goodness of the reproduction of the data, paying no 
attention to the characteristics of the effect in terms of the depth z, makes little sense. 

The measurements in the area corresponding to the penumbra are very sensitive to 
the (input) value used for the lateral displacement of the block; small inaccuracy in this 
value affects the description of the data significantly, introducing spurious effects in the 
function. This area, albeit very important in the determination of the value of the 
parameter A, was not included in this part of the analysis. 

The resulting values are given in Table [3l separately for the four z positions at 
which the measurements have been obtained; it is evident that, for all depth values, the 
quality of the reproduction of the measurements when including the block-scattering 
effects is superior to the case that only the pristine-beam contribution is considered. 
The importance of the inclusion of these effects decreases with increasing distance to 
the block. The improvement for z = 150 mm when including the block-scattering effects 
is impressive. Judging from the values for the given degrees of freedom, there can be 
no doubt that the overall reproduction of the data is satisfactory. Last but not least, 
despite the fact that the description of the option-8 measurements is debatable, our 
conclusions do not depend on the treatment of the data in that option (i.e., inclusion 
in or exclusion from the database). 

Out of the available 120 lateral fluence profiles, corresponding to z = 150 mm, only 
five profiles did not show improvement after the scattering effects were included. In two 
of these profiles, namely the energy-NeT combinations of (176.08 MeV, 102.2 mm) and 
(179.00 MeV, 103.3 mm), the scattering contributions are present in the measurements, 
yet at different amounts compared to the MC-generated data; additionally, a hard-to- 
explain slope (i.e., of different sign to what is expected after including the scattering 
contributions) is clearly seen in these measurements around y = mm. On the 
other hand, no scattering effects ('ears') are discernible in the (167.42 MeV, 160.5 
mm), (206.77 MeV, 75.2 mm), and (204.47 MeV, 203.2 mm) combinations at z = 150 
mm. It has to be stressed that these five data sets are surrounded by a multitude of 
measurements showing an impressive agreement between the experimentally-obtained 
and the MC-generated total-fluence distributions; due to this reason, we rather consider 
the absence of improvement (that is, after the scattering contributions are included) in 
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the description of the data in these five profiles as indicative of experimental problems. 

An alternative way of displaying the content of Figs. [20] and [21] is given in Fig. [22] 
instead of showing separately the measurements and the MC data, shown in Fig. [22] are 
the normalised residuals, defined as 



where v^^^ denotes the i-th measurement and vf^^ the corresponding MC-obtained 
fluence; Sv'^^^ and 6vf^'" represent their uncertainties. The advantage of such a plot 
is evident as direct information on the reproduction may be obtained faster than from 
Figs. [20] and [211 fo^' instance, not only can one immediately become aware of the failure 
of the pristine-beam data close to the borders of the block, but also of the severity of 
this failure. Evidently, the pristine-beam contribution underestimates the fluence by 
about 1 to 3 standard deviations for —50 mm < y < —40 mm. It is also interesting to 
note that, after the scattering effects have been included, the normalised residuals Zi 
show significantly smaller dependence on the lateral distance y (ideally, no dependence 
should be seen). 

3. 3. An example of the application of the corrections 

The aim of the present paper was to provide the systematic description of a method 
to be used in the determination and application of the corrections which are due to 
the presence of BL/BSDs in proton planning. Despite the fact that no emphasis was 
meant to be put on a clinical investigation, one simple example (of the application of 
the corrections) may nevertheless be called for. 

To this end, the current version of the Treatment Planning System Eclipse'-"'^^-* 
(Varian Medical Systems Inc., Palo Alto, California) was extensively modified to include 
the derivation (in beam configuration) and the application (in planning) of both block- 
relating corrections; the application of each of the two corrections may be requested 
separately in the user interface. (A review article, providing the details of the dose 
evaluation in proton therapy, as well as its implementation in Eclipse, will appear soon, 
see Ulmer et al (2009).) A simple rectangular water phantom was created, within which 
a planning treatment volume (PTV) of 90 cm^ was arbitrarily outlined. Four one-field 
treatment plans were subsequently created as follows: 

a) a plan without block-relating corrections, 

b) a plan with block-thickness corrections, 

c) a plan with block-scattering corrections, and 

d) a plan with both block-relating corrections % 

jj To void double counting in the case that both block-relating corrections are requested, the block- 
thickness corrections are first applied to the pristine-beam fluence. The block-scattering corrections, 
which already contain their corresponding block-thickness effects, are subsequently added to the 
'thickness-corrected' pristine-beam fluence. 
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In each of these plans, a brass block was inserted and fitted to the cross section of 
the PTV. Subsequently, a dose of 100 Gy was delivered to the PTV, using the double- 
scattering technique of the NCC machine. The block-relating corrections were estimated 
on the basis of = 64; in fact, the results are practically insensitive to the value of 
A'^, for > 16. Finally, the resulting dose maps were compared; the dose differences 
of plans (b), (c), and (d) to plan (a) were estimated and compared (Figs. [23] and [2^ . 
leading to the following conclusions. 

• As expected, the application of the block-thickness corrections results in lower dose 
values. This is due to the fact that part of the incident flux is blocked as a result 
of the nonzero thickness of the block. 

• As expected, the application of the block-scattering corrections leads to higher dose 
values. This is because some protons, which would otherwise fail to contribute to 
the fluence (as they impinge upon the block), scatter off the material of the block 
and 're-emerge' at positions in the bore or on the downstream face of the block. 

• As far as the delivered dose is concerned, the effects of the block thickness and 
block scattering 'compete' one another. In the water phantom used in this section, 
the block-thickness corrections dominate. 

• The presence of the block in the plan of the water phantom used in this section 
induces effects which amount to a few percent of the prescribed dose (see Fig. IMl) . 
The largest effects appear in the area neighbouring the border of the block. It 
is also worth noticing the characteristic contributions of the scattered protons in 
the entrance region in the frontal and sagittal views in Fig. [2H the largest part of 
the dose in the entrance region corresponds to the low-energy component of the 
scattered beam. As the local dose, delivered in the entrance region, is significantly 
smaller than the corresponding value delivered to the target (which, in fact, is an 
agrument in favour of the use of protons in radiation therapy), the corrections which 
one has to apply to it, though representing a small fraction of the prescribed dose, 
are sizable. 

4. Conclusions 

The present work deals with corrections which are due to the presence of beam-limiting 
and beam-shaping devices in proton planning. The application of these corrections 
is greatly facilitated by decomposing the effects of two-dimensional objects into one- 
dimensional, easily- calculable contributions (miniblocks). 

In the derivation of the thickness corrections, we follow the strategy of Slopsema 
and Kooy (2006). Given the time restrictions during the planning, the derivation of the 
scattering corrections necessitates the introduction of a two-step approach. The first 
step occurs at the beam-configuration phase. At first, the value of the only parameter 
of our model (A) is extracted from the half-block fluence measurements. A number 
of Monte-Carlo runs follow, the output of which consists of the parameters pertaining 
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to convenient parameterisations of the fluence contributions of the scattered protons. 
These runs take account of the variabihty in the block material and thickness, incident 
energy, and nozzle-equivalent thickness in all the options for which a proton-treatment 
machine is configured. To enable the easy use of the MC results, the output is put in 
the form of expansion parameters in two geometrical quantities which are involved in 
the description of the scattering effects. The scattering corrections for all the blocks in a 
particular plan are determined from the results, obtained at beam-configuration phase, 
via simple interpolations. 

The verification of the method should involve the reproduction of dedicated dose 
measurements. At present, given the lack of such measurements, the only possibility for 
verification rested on re-using the half-block fluence measurements, formerly analysed 
to extract the A value; this is a valid option because parts of the input data had 
been removed from the database to suppress the (present in the measurements) block- 
scattering contributions. We investigated the goodness of the reproduction of the 
measurements on the basis of the function and concluded that the inclusion of the 
scattering effects leads to substantial improvement. 

The method presented in this paper was applied to one plan involving a simple 
water phantom; the different contributions from the two block-relating effects have 
been separately presented and compared. These effects amount to a few percent of the 
prescribed dose and are significant in the entrance region and in the area neighbouring 
the border of the block. 
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Table 1. The weighted averages of the extracted A values (and their statistical 
uncertainties) for the eight options of the NCC machine. The uncertainties are shown 
only for the sake of completeness; they have not been taken into account in the results 
of Section [3. 1.31 



Option number 


A(5A) 


1 


4.667(0.046) 


2 


5.020 (0.070) 


3 


4.549 (0.057) 


4 


4.857(0.069) 


5 


4.800 (0.090) 


6 


5.26 (0.11) 


7 


6.42 (0.15) 


8 


4.65(0.27) 



Table 2. The absolute yields (numbers of particles) of the different types of the 
scattered protons for an aperture size of 100 mm, for three energy-NeT combinations 
of the NCC machine (see text); a 65-mm thick brass block has been used. 



E^ax (MeV) 


OTs 


BSITs 


GTITs 


79.5 


35696 


64130 


23799 


150.2 


72015 


42128 


54383 


212.8 


319806 


83808 


242229 



Table 3. The values, corresponding to the reproduction of the half-block fluence 
measurements, separately for the four z positions at which the data have been obtained; 
evidently, the quality of the reproduction of the measurements when including the 
block-scattering effects is superior to the case that only the pristine-beam contribution 
is taken into account. NDF denotes the number of degrees of freedom. The lower part 
of the table contains the results after excluding the measurements of option 8. 



z (mm) 


Pristine 


Total 


NDF 


150 


22710.25 


8551.34 


14511 





15384.09 


11826.71 


14210 


-150 


13229.26 


12255.67 


13841 


-300 


12949.22 


12573.03 


13628 


150 


16024.44 


5761.21 


12704 





12001.81 


9761.32 


12404 


-150 


10751.43 


10277.58 


12057 


-300 


10658.78 


10531.18 


11843 
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Figure 1. The outer tracks (OTs) and the bore-scattered inner tracks (BSITs) emerge 
from the bore. The going-though inner tracks (GTITs) emerge from the downstream 
face of the block. 




Figure 3. The effects of the BL/BSD have to be evaluated at the point Q, which is 
projected to the point P on the BL/BSD plane. In this figure, the point P lies within 
the area B. 
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Figure 4. Evaluation of tiie effects of tlie BL/BSD at the point Q of Fig. [3] (not 
shown), whose projection onto the BL/BSD plane is the point P, on the basis of four 
directions (resulting, in this case, in four miniblocks) . In this figure, the point P lies 
within the area B. 



Collimator effects in proton planning 



30 




Figure 5. Evaluation of tiie effects of tlie BL/BSD at the point Q of Fig. [3] (not 
shown), whose projection onto the BL/BSD plane is the point P, on the basis of four 
directions (resulting, in this case, in five miniblocks). In this figure, the point P lies 
outside the area B. 
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Figure 6. The relation between the DICOM block, the extension, and the downstream 
projection on the (x,y) plane (calculation plane) at a specified depth z. 




extension 

downstream projection 

upstream projection 



Figure 7. The extension and the projections of the downstream and upstream faces 
of the block on a calculation plane. The point P lies within the extension, the point 
P' without. 




Figure 8. The essentials for the evaluation of the contribution of a miniblock to the 
fluence at specified points (e.g., at P and P') on the calculation plane. 
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Figure 9. Due to onc-diniensional nature of the miniblocks, a point lying within the 
extension (of a miniblock) may also receive contributions which are characteristic to 
points lying in the exterior of the extension. 




Figure 10. The description of the kinematics inside the block. The auxiliary 
coordinate system introduced in this figure should not be confused with the formal 
coordinate system of Fig. [S] 
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Lateral fluence distributions 
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y/mm 



Figure 11. The lateral fluence distributions for one energy-NeT combination (171.99 
MeV, 154.5 mm) of option 2 of the NCC machine. 
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Figure 12. The values of the parameter A for all the energy-NeT combinatfons of all 
double-scattering options of the NCC machine. 
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Figure 13. A typical description of the lateral fluence distribution for outer tracks. 
Bins with fewer than 10 entries are not shown. 
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Inner tracks (BSITs and GTITs) 



• MC-generated data o Fitted data 
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Figure 14. A typical description of the lateral fluence distributions for inner tracks. 
Bins with fewer than 10 entries are not shown. 
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z=100 mm 
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Figure 15. One case of the scattering corrections (to be applied to the lateral fluence 
distribution of the pristine tracks) at ^; = 100 mm. The lateral displacement of the 
block was 6 = 20 mm. 
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z=0 mm 
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Figure 16. One case of the scattering corrections (to be applied to the lateral fluence 
distribution of the pristine tracks) ai z = mm (i.e., at isocentre). The lateral 
displacement of the block was b = 20 mm. 
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z=-100 mm 



AOTs OGTITs •BSITs 




-30 -20 -10 10 20 30 40 50 60 70 



y/mm 



Figure 17. One case of the scattering corrections (to be applied to the lateral fluence 
distribution of the pristine tracks) at ^ = —100 mm. The lateral displacement of the 
block was b = 20 mm. 
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Energy distribution of the outer tracks 
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Figure 18. The energy distributions for the outer tracks for three energy-NeT 
combinations of the NCC machine (see text). 
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Energy distribution of the inner tracks 
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Figure 19. The energy distributions for the inner tracks for three energy-NeT 
combinations of the NCC machine (see text). 
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E=1 58.42 MeV, NeT=125.3 mm, z=150 mm 



o Monte-Carlo data: pristine tracl<s Half-blocl< measurements 
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Figure 20. The lateral fluence measurements (continuous line) corresponding to one 

energy-NeT combination of one option of the NCC machine, taken 100 mm away from 
the downstream face of the block. The Monte-Carlo data shown correspond only to 
the pristine-beam fluence obtained at the same incident-energy, NeT, and z values; 
the measurements have been scaled up by a factor which is equal to the ratio of the 
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E=1 58.42 MeV, NeT=125.3 mm, z=150 mm 

o Monte-Carlo data: all tracks ^—Half-block measurements 
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Figure 21. The lateral fluence measurements (continuous line) corresponding to one 
cncrgy-NeT combination of one option of the NCC machine, taken 100 mm away from 
the downstream face of the block. The Monte-Carlo data shown correspond to the total 
(pristine-beam plus scattered-protons) fluence obtained at the same incident-energy, 
NeT, and z values; the measurements have been scaled up by a factor which is equal 
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E=1 58.42 MeV, NeT=125.3 mm, z=150 mm 



o Pristine tracks • All tracks 



(0 
3 

■a 

'in 

Oi 

■a 

E 
o 



o 



o o 
6> ° 



• o • 



o ^» 

o 



B 8 gO^r^ ^ .O « ® 'P-^^ O 



o 
o 



•o 



" CP 

o o« 

o ^ 

-fP— ' — OH 



o o 



. «8 



o o 

|o 



o cp 



-50 



-40 



-30 



-20 

y/mm 



-10 



10 



Figure 22. An alternative way of displaying the contents of Figs. [^D] and [^T] shown in 
this figure are the normalised residuals, plotted versus the lateral distance y. Evidently, 
the pristine-beam contribution underestimates the fluence by about 1 to 3 standard 
deviations for —50 mm < y < —40 mm. After the inclusion of the block-scattering 
effects, the residuals nicely cluster around 0. 
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Figure 23. The dose contributions corresponding to the block-thickness (on the left) 
and to the block-scattering (on the right) corrections for the simple water phantom of 
Section [531 
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Figure 24. The dose contributions corresponding to both block-relating corrections 
for the simple water phantom of Section [3.31 



